#!/bin/awk -f
BEGIN{
# bond type 1: tiling edges
  bd1[0]="10000"; 
  bd1[1]="01000"; 
  bd1[2]="00100"; 
  bd1[3]="00010"; 
  bd1[4]="00001";
  bd1[5]="-10000"; 
  bd1[6]="0-1000"; 
  bd1[7]="00-100"; 
  bd1[8]="000-10"; 
  bd1[9]="0000-1";
# bond type 2: perpendicular bonds
  bd2[0]="1-1000"; 
  bd2[1]="01-100"; 
  bd2[2]="001-10"; 
  bd2[3]="0001-1"; 
  bd2[4]="-10001";
  bd2[5]="-11000"; 
  bd2[6]="0-1100"; 
  bd2[7]="00-110"; 
  bd2[8]="000-11"; 
  bd2[9]="1000-1";
# hamiltonian
  j1=0.01; j2=0.02;
# temperatures
  beta="0.1 3 6 8 10 12 14 16 18 20 50 100";
  nloops=1000; emin=1000000;
  for(i in ARGV){
    if(ARGV[i]=="-h"){
      print "Synopsis : mctspin [conf=011100...] nl=LOOPS beta=b1,b2,b3... j1=J1 j2=J2 TILINGNAME";
      print "** Default init conf=000000... ";
      print "** Default schedule beta=("beta")";
      print "** Default loops foreach beta nl="nloops" (is mult by Nnodes)";
      hexit=1; exit;
    }
    m=split(ARGV[i],a,"=");
    if(m==2){
      if(a[1]=="nl"){nloops=a[2]};
      if(a[1]=="beta"){gsub(/,/," ",a[2]);beta=a[2]};
      if(a[1]=="conf"){split(a[2],conf,"")};
      if(a[1]=="j1"){j1=a[2]};
      if(a[1]=="j2"){j2=a[2]};
    }else{
      fn=ARGV[i];
    }
  }
  ARGC=2; ARGV[1]=fn;
}
FNR==1{b5a=$1" "$2" "$3" "$4" "$5;split(b5a,ba)} # basis vector 1 (5D)
FNR==2{b5b=$1" "$2" "$3" "$4" "$5;split(b5b,bb)} # vasis vector 2 (5D)
FNR==3{ # number of tiling nodes
  mu=$2;
  if($2>1){print "This is stacked tiling.";hexit=1;exit}
  nn=$1;
  if(conf[1]==""&&$3~/spinconf/){
    split($3,a,"=");
    split(a[2],conf,"");
  }
}
FNR>3{
  n++;
  t[n]=$1" "$2" "$3" "$4" "$5;
  if(conf[n]==""){s[n]=0}else{s[n]=conf[n]}
}
END{
  if(hexit==1){exit};
  print "tiling="fn" J1="j1" J2="j2;
  eb["100"]=j1;
  eb["110"]=-j1;
  eb["101"]=-j1;
  eb["111"]=j1;
  eb["200"]=j2;
  eb["210"]=-j2;
  eb["201"]=-j2;
  eb["211"]=j2;
## find bonds, looping over nodes AND cell images
  for(i=1;i<=n;i++){ # node loop 1
    split(t[i],p);
    for(j=i+1;j<=n;j++){ # node loop 2
      split(t[j],q);
      for(k=-1;k<=1;k++){ # cell image loop 1
	for(l=-1;l<=1;l++){ # cell image loop 2
	  for(m=1;m<=5;m++){
	    b[m]=q[m]-p[m]+k*ba[m]+l*bb[m];if(b[m]=="-0"){b[m]=0}
	  }
	  ss=b[1]b[2]b[3]b[4]b[5]; # plain string representation of the bond
	  for(d in bd1){
	    if(ss == bd1[d]){
	      b[i" "j" "ss]=1; b[j" "i" "ss]=1;
	      bl[i]=bl[i]" "j":1";
	      bl[j]=bl[j]" "i":1";
	    }
	  }
	  for(d in bd2){
	    if(ss == bd2[d]){
	      b[i" "j" "ss]=1; b[j" "i" "ss]=1;
	      bl[i]=bl[i]" "j":2";
	      bl[j]=bl[j]" "i":2";
	    }
	  }
	}
      }
    }
  }
  e=0; for(i in t){e+=energy(i)};
  e=e/2; # we have double-counted all bonds
  nt=split(beta,be);
######### Could implement sometime if needed to speed up
# initialize table of exponentials [inside each beta loop]
#  bmx=6;
#  for(i=-2*nbmx;i++;i<=2*nbmx){
#    for(j=-2*nbmx;j++;j<=2*nbmx){
#      arg=i*j1+j*j2;
#      if(arg>=0){et[arg]=exp(-be[x]*arg)}
#    }
#  }
###########################################################
  printf "init conf="; for(i=1;i<=nn;i++){printf s[i]}; print " E="e;
  for(it=1;it<=nt;it++){
    na=0;
    for(l=1;l<=nn*nloops;l++){
      ii=int(nn*rand())+1;
      eo=energy(ii); # energy for THIS node
      s[ii]=1-s[ii]; # swap state of the node
      en=energy(ii); # new energy
      de=en-eo; # note: we dont doublecount here!
#      print ii,s[ii],de;
      if(de<=0){ # accept
	e=e-eo+en; na++;
      }else if(de>0){
	x=rand();
	earg=de*be[it];
	if(earg>99){expo=0}else{expo=exp(-earg)};
	if(expo>x){ # accept+
	  e=e-eo+en; na++
	}else{ # reject : swap state backward
	  s[ii]=1-s[ii];
	}
      }
      e=1.0*sprintf("%1.10f",e);
      if(e<=emin){
	if(e<emin){split("",sco)};
	emin=e;
	cf=""; for(i=1;i<=nn;i++){cf=cf""s[i]};
	sco[cf]++;
      };
    }
    ar=na/nn/nloops;
    printf("%2i %5.1f %6.2f %8.6f\n",it,be[it],e,ar);
  }
  printf "final conf="; for(i=1;i<=nn;i++){printf s[i]};
  en=0; for(i in t){en+=energy(i)}; en=en/2;
  print " E="e" E_check="en;
  sub(/^.*\//,"",fn);
  print "Lowest-E ensemble for E_min="emin" ";
## eliminate inverse states
  j=0;for(i in sco){
    sinv=i;
    gsub(/1/,"2",sinv);
    gsub(/0/,"1",sinv);
    gsub(/2/,"0",sinv);
    if(sce[i]==""&&sce[sinv]==""){sce[i]++}; # counts are not correct now
  }
  j=0;for(i in sce)if(sce[i]!=""){
    j++; 
    fo=fn".cfg"j;
    print fo,i;
    print b5a>fo;
    print b5b>>fo;
    print nn" "mu" # spinconf="i>>fo;
    split(i,a,"");
    for(k=1;k<=nn;k++){print t[k],"  "a[k]>>fo}
  }
}

function energy(ii,  b,c,n,i,e,x){
  n=split(bl[ii],b);
  for(i in b){
    split(b[i],c,":");
    x=eb[c[2]s[ii]s[c[1]]]; 
    if(x==""){print "energy-tab-error"}else{e+=x};
  }
  return e;
}
